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Planar and curved droplet networks 
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We study how a two-dimensional square sheet of droplets evolves to the preferred triangular lattice. 
The rearrangement can be induced either by an impurity seed for flat sheets or by curving the sheet. 
We analyze the propagation of transformation fronts and the spatial reordering subsequent to the 
initial perturbation. 
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I. INTRODUCTION 

Metastable states, local minima of the free energy land¬ 
scape, are very common. Transitions from one local 
minimum to another local minimum (not necessarily the 
ground state) may take place spontaneously, or when sys¬ 
tem parameters, such as the ambient temperature and 
pressure, are changed. They may also be triggered by 
external perturbations or the introduction of impurities. 
Examples include avalanches (T|, protein folding [2], de¬ 
vitrification of glass [3j, crystallization of supercooled liq¬ 
uids [1| and aggregation and sedimentation of colloidal 
systems [5j. One broad class of such transitions is solid 
to solid phase transitions. Diamond, for example, is not 
thermodynamically stable at standard temperature and 
pressure - it exists because of the large activation bar¬ 
rier for conversion to graphite [6J. Martensite, which re¬ 
sults from rapid cooling, is also a metastable phase and 
is easily destroyed by heating [7j. Despite their physical 
significance the detailed dynamics of solid to solid phase 
transitions is difficult to observe experimentally as it is 
rapid and typically occurs on small length scales. Col¬ 
loidal crystals, whose dynamics are significantly slower, 
provide a rich analogue system to their atomic or molec¬ 
ular counterparts mm- 

An entirely different model system is that of droplet 
networks HU formed by printing tens of thousands of 
micron-sized aqueous droplets that bind by forming sin¬ 
gle lipid bilayers mm to form a cohesive network in 
bulk oil. Functionality can be added to the network and 
the droplets can be programmed with different osmolari- 
ties (osmolarity is the measure of solute concentration) to 
fold the network into certain designed structures hu mi. 
Here instead of discussing the promising applications of 
the network, we will explore the potential possibilities of 
using the droplet network as a platform to study lattice 
transitions. Both the colloidal crystals and the droplet 
systems are platforms to study lattice transitions: Col¬ 
loidal crystals have long been studied, while the droplet 
systems have their own advantages such as well-defined 
interaction HU, programmable initial states HU and re¬ 
versibility of a folded shape [14] • As the initial posi- 
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tions of the droplets are predetermined in this system, 
the droplets can be arranged in metastable states. This 
is an ideal system to study lattice transitions. The in¬ 
teraction between droplets may be simply modeled as a 
spring potential with a damping term HU- We exam¬ 
ine a simple metastable state - a two-dimensional (2D) 
square network - and study its transformation into a tri¬ 
angular lattice. This transformation has also been inves¬ 
tigated in a study of temperature-induced lattice transi¬ 
tions in colloidal crystals [15] . Osmolarity gradients lead 
to flow of water from lower to higher osmolarity, result¬ 
ing in the swelling of the higher osmolarity droplets and 
the shrinking of the lower osmolarity droplets. We thus 
use osmolarity gradients to introduce size polydispersity. 
The effects of size polydispersity have been analyzed for 
a variety of physical processes, such as crystallization 
m and granular dynamics HU- Yao et al. m sug¬ 
gest that polydispersity-driven topological defects could 
act as order restoring excitations. We start by consid¬ 
ering the situation where there is an impurity of larger 
size in the square lattice. Surrounded by smaller par¬ 
ticles, the coordination number of the impurity tends to 
exceed four, destroying the regular square lattice and ini¬ 
tiating the rearrangement. We focus on the dynamics of 
the transformation, rather than minimal energy states 
HMU, with the final configurations not necessarily the 
ground states. After the rearrangement is initiated, we 
explore how the subsequent propagation of reordering re¬ 
lates to the strengths of the interaction and the damping. 
For two initial impurities we find that the competition 
between transformation fronts, as the reordering propa¬ 
gates, leads to the formation of an intervening domain 
wall. 

We also consider the motion of droplets on curved sur¬ 
faces, motivated by previous studies of colloidal parti¬ 
cles confined on a two-sphere (22 as well as on capillary 
bridges [23l - f25j . and explore curvature-induced transfor¬ 
mations of the lattice configuration of the droplet net¬ 
work. This may be realized experimentally by choosing 
appropriate hydrophilic substrates. In the quasi-static 
process of slowly increasing the curvature of an elliptic 
paraboloidal surface we observe a first transition after 
which only the boundary regions transform into the trian¬ 
gular lattice structure and a second transition after which 
a cross shape forms in the central area. Rearrangements 
take place in stages also for the case of a paraboloidal 
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surface. For impurity-driven transformations on a spiral 
surface we find that ordered domains are rotated in the 
final configuration. 


II. MODEL AND METHODS 

We briefly describe the model of Ref. m- Lipid-coated 
aqueous droplets in a continuous oil phase adhere at their 
interfaces by forming stable bilayers. To simulate this 
mechanical interaction, every droplet i is treated as a 
point mass with an associated radius Ri. If two droplets 
i and j come within a distance Ri + Rj, they become 
connected by a Hookean spring of equilibrium length 
a(Ri + Rj): the spring represents the attractive interac¬ 
tion between fused droplets and the parameter a, which 
is less than 1 , approximates droplet deformation due to 
droplet-droplet adhesion. The net force on each droplet 
i is therefore 


Fi= ^2 k [ rij - a (Ri + Rj)] fij - 7 $, ( 1 ) 

oeCi 

where C,; is the set of droplets connected to the ith 
droplet, fij is the vector from droplet i to droplet j with 
magnitude r t j , k is a spring constant, 7 is a damping 
coefficient, and v, is the velocity of the zth droplet. The 
force increases linearly with the distance between two 
droplets until the distance reaches Ri + Rj, where the 
force drops to zero and remains zero for distance above. 
The position of each droplet is updated using 


Fi = m 


dt 2 


( 2 ) 


where m is the mass of the droplet. The equation of 
motion for a droplet is given by 


d 2 fj 

dt 2 


T, w 2 [r^ - a (Ri + Rj)] fj - u>(fi, 

jeCi 


(3) 


where u> = \J~kJm is the angular frequency and £ = 
7 /y/mk is the dimensionless damping coefficient. 


III. RESULTS AND DISCUSSION 
A. Impurity-induced transformation 

Defects play an essential role in understanding 2D crys¬ 
talline systems. Defects on a shell, for example, may 
drive the shape from round to faceted, depending on the 
Foppl-von Karrnan number [26], and the precise defect 
distribution affects the shell’s ability to sustain external 
pressure m- Here we use osmolarity gradients to intro¬ 


duce impurities, and analyze how the subsequent prop¬ 
agation of rearrangements depends on the values of the 
spring constant k and the damping parameter £. Osmo¬ 
sis is implemented as follows: when two droplets with 
different osmolarities are connected by a spring, the vol¬ 
ume of water transferred per unit time from a droplet i 
with osmolarity C, to a droplet j with osmolarity Cj is 
given by Fick’s first law as 


Jij = AijD (Cj — Ci) , (4) 

where A ,; 7 is the area of the bilayer between the two 
droplets, and D denotes a permeability coefficient which 
is assumed to be constant and identical for all bilayers. 
Following Ref. E5. we assume k, m and D are the same 
for all droplets. 

The positions of the droplets are updated in the fol¬ 
lowing way: at time t, any two droplets i and j that 
have come within a distance Ri + Rj are connected by a 
spring of natural length a(Ri + Rj). If the osmotic inter¬ 
action is active, the volume of water transferred between 
each pair of droplets joined by a bilayer is calculated 
according to Eq. (4), and the size of each droplet is up¬ 
dated accordingly. The position of each droplet at time 
t + At is calculated according to Eqs. (1) and (2) using a 
fourth-order Runge-Kutta scheme. All the droplets have 
uniform initial volume Vq = 47 t.Rq/ 3 , where Ro is the 
initial radius. We set R 0 = 1, m = 1, a = 0.8, D = 1 
(D is chosen larger than that in Ref. [TT| to shorten the 
time required for the initial water exchange - this has 
little effect on the subsequent propagation of rearrange¬ 
ments), At = 10~ 3 (in the large damping case Fig. If, 
At = 2 x 10~ 4 to make the iteration stable). The units 
of length, mass, spring constant and osmolarity can be 
obtained by matching their experimental values given in 
Ref. [11]; the units of other quantities are a combination 
of these units. 

We start by studying the situation of one impurity. 
Figure 0 ' shows the initial configuration of N = 1225 
droplets positioned in a square lattice, with one impu¬ 
rity of high osmolarity 5/Vo at the center and all oth¬ 
ers of uniform low osmolarity 1/Vo. These droplets are 
first equilibrated without water exchange. To minimize 
boundary effects we implement absorptive boundary con¬ 
ditions by gradually increasing the £ values of the out¬ 
ermost six layers. The setup is illustrated in Fig. [I]): 
the internal droplets have £ = 1 the six outermost lay¬ 
ers have £ = 2,4,10,10,20,20 in order of innermost to 
outermost. With these £ values we do not observe any 
reflective propagations. The boundary conditions when 
the internal droplets have other £ values are given in 
Appendix A. To better describe the system, we define 
the connectivity of the system as the structure of the 
spring network at a given time. After water exchange 
is turned on the impurity swells and connects with the 
four next-to-nearest neighbors - the coordination num¬ 
ber of the impurity seed thus doubles from four to eight. 
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The four next-to-nearest neighbors acquire coordination 
number five. This initiates the propagation of a cross¬ 
shaped region of triangular lattice within the original 
square lattice. Away from the impurity, the additional 
kinetic energy comes from the potential energy of the ini¬ 
tial system through the formation of extra bonds and the 
propagation speed depends on the stiffness (k) and the 
damping (£). Either higher stiffness or lower damping £ 
increase the conversion rate to the triangular lattice. We 
plot the total number of edges of the converted triangu¬ 
lar lattice, for a variety of stiffnesses and damping values, 
versus time in Fig. [1} a nd the corresponding connectivi¬ 
ties at t = 20 in Fig. Ill (c-h). In our simulation time the 
rearrangement is localized in the case of small stiffness 
and high damping (Fig. [TJ: and movie 00 - For higher 
stiffness or lower damping the coarsening transformation 
front propagates (Fig. [l] d-h and movie [TJ d-h). The 
propagation speed is almost constant, as can be seen in 
Fig- 0 , until the front reaches the boundary region where 
it slows down. Figure [T]l will propagate to the bound¬ 
ary provided longer simulation time. Triangular domains 
and square domains are separated by lines of fivefold de¬ 
fects. The position of the defect lines may depend on k , 
£ and the initial perturbation. 

We note that both k and £ may depend on the temper¬ 
ature, and the effect of temperature can be included as 
we can change k and £ for different temperatures. In the 
following we choose k = 10 3 and £ = 0.2, similar to the 
values used in Ref. m- We adopt absorptive boundary 
conditions for impurity-induced transformations. 

We treat now the case of two impurities. Figure 
shows the initial configuration - a square lattice of N = 
2905 droplets with a pair of impurities of high osmolarity. 
The distance between the pair of impurities is 42 lattice 
spacings. Figure [2 ]d shows the connectivity of this system 
at time t = 10. Coarsening fronts continue to propagate 
to the boundary after meeting in the region between the 
impurities, leading to the formation of a domain wall with 
vacancies (Fig. |2j c and d; movie [2j c and d). Different 
final states could be obtained if we arrange the positions 
of the two impurities differently or if we consider more 
impurities. 




B. Effects of randomness 

As there may be stochastic noise in the experiment, 
we study the effects of the randomness of the droplet 
positions. In Fig. |3j we give the droplets in the cen¬ 
tral region (inside the yellow dashed line) a perturba¬ 
tion in a random direction, with a maximal value d, as 
the initial state. Then we relax the system. When d 
is small, the connectivity of the square lattice is pre¬ 
served. After relaxation, the droplets tend to return to 
their equilibrium positions, as shown in Fig. [3^(right), 
where we see a nearly perfect square lattice connectivity. 
Note that the final state may be slightly different from a 
perfect square lattice even after running the simulation 


FIG. 1. (Color online) (a) Initial osmoloarity profile of a 
square lattice in a square area, composed of N = 1225 = 
(Nl + Nb) x ( Nw+Nb ) = (23 + 12) x (23 + 12) droplets, with 
Nl, Nw, Nb being the number of droplets along the length, 
along the width, in the boundary (outside the black dashed 
line), respectively. An impurity of high osmolarity 5/Vo (red) 
at center and others of uniform low osmolarity 1/Vo (blue), 
(b) The absorptive boundary condition in the case £ = 1. 
Color represents £ value. Internal droplets are with £ = 1. 
From innermost to outermost, £ = 2,4,10,10,20,20, respec¬ 
tively, for the six layers in the boundary, (c-h) Connectivity 
at t = 20. Fivefold vertices are shown in red, sixfold in blue, 
sevenfold in yellow and eightfold in grey, (c-e) Fixed k = 1, 
£ = 10,1,0.1. (f-h) Fixed k = 10 3 , £ = 10,1,0.1. (i) Num¬ 
ber of edges with triangular order for various k and £ values. 
Snapshots were generated using the Visual Molecular Dynam¬ 
ics (VMD) package [28] and rendered using the Tachyon ray 
tracer (29] . 

































































































































































































































































4 



FIG. 2. (Color online) (a) Initial configuration of IV = 2905 
(Nl = 83, Nw = 35, boundary not included) droplets with 
two impurities of high osmolarity. The number of columns 
between the pair of impurities is n p = 41. Colors represent 
osmolarity, as in Fig. la. (b) Connectivity at t = 1.4. (c) 
Connectivity at t = 10. (d) Connectivity at t = 10 for a 

Nl = 82, Nw = 35, n p = 40 system. 


for a rather long time (say t ss 50T, with T the period 
and given by T = 2tt/uj) due to the damping term. On 
the contrary, when d is large enough, the initial square 
lattice will break down without introducing additional 
defects. After relaxation, the transformation from the 
square lattice to the triangular lattice will spread out 
of the core region and reach the boundary, as shown in 
Fig. [3jo(right) . That the effects of randomness depend 
on strength is understandable as the square lattice is a 
metastable state and thus has the ability to sustain small 
perturbations. Provided the noise in an experimental set¬ 
ting is not large and the initial positions of the droplets 
are carefully arranged, and subsequently allowed to re¬ 
lax for a sufficiently long time, we expect our simulation 
results to provide an accurate guide. 




FIG. 3. (Color online) Left: Initial configurations with small 
perturbation d = 0.13 (a) and larger perturbation d = 0.3 
(b) in the central regions (the yellow dashed line). The 
perturbation can be seen from the irregular white space be¬ 
tween droplets. The black dashed line indicates the boundary. 
Right: Connectivity at t = 10 after relaxation. 


C. Curvature-induced transformation 


Previous studies have explored chemically-driven free- 
running droplets on solid substrates [50H55] . Here we 
shift to treating the case of droplets constrained to move 
on curved surfaces which can arise due to osmolarity- 
gradient induced bending of an initially flat surface. For 
simplicity we assume that droplets on the surface are not 
subject to any pressure (we don’t consider the motion 
in the normal direction). With this assumption, it can 
be shown that the the differential equations of motion of 
a droplet of mass to, constrained to move in a smooth 
surface S: x(u, v), read (T?X 

u = (j) — r^w 2 — 2t\ 2 uv — (5) 

and 

v = ip- - 2Tj 2 uv - r 2 2 u 2 , (6) 

where P^ (a,/3,7 = 1,2) are Christoffel symbols [35], (p 

and ip can be written in term of the force vector F de¬ 
fined in the Euclidean space and given by Eq. (1) as <p = 
[ 922 (F ■ dx/du )^- gi 2 {F ■ dx/dv)\/mg and ip = bffi 2 (-F • 
dx/du) + gn(F ■ dx/dv)]/mg , with g a p (a, (3 = 1,2) the 
metric tensor and g = det (g a p) = gng 2 2 — (gi 2 ) 2 - The 
interaction between droplets and the damping effect are 
included in the <p and ip terms; the curvature effect is 
included in the other three terms which vanish in flat 
geometry in each equation. 
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FIG. 4. (Color online) N = 1245 droplets constrained to 
move on an elliptic paraboloid surface (a) and a hyperbolic 
paraboloid surface (b). (a) Connectivities after the transitions 
at a ~ 8.9 (left) and at a ~ 7.2 (right), (b) Connectivities 
after the transitions at a ~ 8.9 (left) and a ~ 8.4 (right). 


We start by studying a curvature-induced transfor¬ 
mation. N = 1245 droplets of uniform osmolarity are 
positioned within a circular area yj x 2 + y 2 < 20 on an 
elliptic-paraboloid surface given b y z = x 2 /a 2 + y 2 /a 2 , 
with a = 15, and subsequently equilibrated. For this 
elliptic-paraboloid surface, F] 1 = V\ 2 = Ax/ (a 4 + Ax 2 + 
Ay 2 ), T 2 n = r| 2 = Ay/(a* + Ax 2 + Ay 2 ), r| 2 = r? 2 = 0. 
To simulate a quasi-static process we relax the system 
for 2 x 10 3 time steps every time a is decreased by 0.1. 
We observe a first transition (movie [4^) at a « 8.9, with 
boundary regions transforming into a triangular lattice 
with defects, while the central region preserves the square 
lattice structure. The connectivity after this transition is 
shown in Fig. [4ji (left). As a further decreases, we find a 
second transition happens at a sa 7.2, with the formation 
of a cross shape in the central region (Fig. EK right). 
After this transition the remaining square lattice grad¬ 
ually transforms into triangular lattice as a decreases. 
The situation is much the same for droplets confined on 
a hyperbolic paraboloid surface, z = x 2 /a 2 — y 2 /a 2 with 
a = 15. We observe two transitions (movie |4 ]d), the first 
at o » 8.9 (Fig. §>> left) and the second at a ~ 8.4 
(Fig#, right). Note that the surfaces chosen here are 
commensurate with the 4-fold symmetry of the original 
square lattice and we don’t include randomness. After 
each transformation the lattice still preserves the 4-fold 
symmetry. Curvature-induced transitions might be seen 
in systems with a substrate of non-zero Poisson ratio, as 
a region of non-zero Gaussian curvature is likely to form 
upon stretching the substrate along the axial direction. 


(a) 




FIG. 5. (Color online) (a) The spiral surface given by Eq. [7]. 
(b) Connectivity (boundary not shown) at t = 10, with the 
initial configuration same as in Fig. la but on the spiral sur¬ 
face. 


D. Impurity-induced lattice transformations on a 
spiral surface 

We next explored the lattice evolution with both im¬ 
purities and spatial curvature. We position the system 
described in Fig. E* on a spiral surface (Fig. io given by 

* = K*) 2 + (-) 2 ] 1 sin[4arctan(-) + 4</(—) 2 + (*#7) 
a a x V a a 

with a = 20, and equilibrate it without exchange of wa¬ 
ter. Figure [5 |d (movie i>) shows the connectivity at 
t = 10 after the exchange of water is turned on. The 
configuration exhibits a rotation of domains compared 
to its flat analogue. This is a simple example showing 
how curvature may affect the propagation of the rear¬ 
rangement induced by one impurity; other patterns may 
be obtained by various curvature designs and impurity 
arrangements. 


IV. CONCLUSION 

We have studied a 2D droplet system with an inter¬ 
action described by a spring potential with a damping 
term. To explore the potential possibility of using the 
droplet network to study lattice transitions, we exam¬ 
ine the transformation from a simple metastable state, 
a square lattice, to a triangular lattice, induced by ei¬ 
ther impurities or spatial curvature. We find that the 
subsequent propagation of rearrangements depends on 
the spring constant and the damping coefficient. When 
randomness is included, it may destroy the square lattice 
even without introducing additional defects if its strength 
is strong enough. Moreover, we show examples that dif¬ 
ferent propagation patterns can be obtained through var¬ 
ious impurity arrangements or curvature structures: the 
competition between transformation fronts induced by 
two impurity seeds leads to the formation of a domain 
wall; the rearrangement takes place gradually from edge 
to central areas if we adapt an elliptic paraboloid or a 
hyperbolic paraboloid surface and increase its curvature 
slowly; an impurity-induced transformation on a spiral 
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surface shows a rotation of domains. 
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V. APPENDIX A. ABSORPTIVE BOUNDARY 
CONDITIONS 

In the case of impurity-induced transformations, we 
add to the system a six-layer boundary with gradually 
increased damping values. For a system with a certain 
damping £, the damping values of the six boundary layers 
(from innermost to outermost) are 2£, 4£, IOC, IOC, 20C, 


20C- Thus for systems with C = 10 (Fig. |T] c and f), the 
C values of the six boundary layers (from innermost to 
outermost) are 20, 40, 100, 100, 200, 200. For systems 
with C = 1 (Fig. [lj d and g), the values are 2, 4, 10, 10, 
20, 20, as plotted in Fig. [TJr. For systems with C = 0.1 
(Fig. [lj e and h), the values are 0.2, 0.4, 1.0, 1.0, 2.0, 
2.0. For systems with ( = 0.2 (Fig. i Fig. [3j Fig. [3] and 
Fig. [5]) , the values are 0.4, 0.8, 2.0, 2.0, 4.0, 4.0. 


VI. APPENDIX B. SUPPLEMENTARY MOVIES 

To better illustrate the transformations, we attach 
Movie [l] (c to h), Movie [2] (c and d), Movie [4] (a and b) 
and Movie [5] which match the number of figures, to the 
article. In the movies, color represent osmolarity, with 
high osmolarity 5/Vo red and low osmolarity 1/Vo blue. 
The gray droplets in Movie [T| represent the absorptive 
boundary. 
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